Phase diagram of the Heisenberg antiferromagnet with four-spin interactions 
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We study the quantum phase diagram of the Heisenberg planar antiferromagnet with a subset of 
four-spin ring exchange interactions, using the recently proposed heirarchical mean-field approach. 
By identifying relevant degrees of freedom, we are able to use a single variational anzatz to map the 
entire phase diagram of the model and uncover the nature of its various phases. It is shown that 
there exists a transition between a Neel state and a quantum paramagnetic phase, characterized by 
broken translational invariance. The non-magnetic phase preserves the lattice rotational symmetry, 
and has a correlated plaquette nature. Our results also suggest that this phase transition can be 
properly described within the Landau paradigm. 

PACS numbers: 75.10.Jm, 64.70.Tg, 75.40.Cx 



I. INTRODUCTION 



Quantum phases of matter and their transitions are 
of fundamental concern to modern condensed matter 
physicsi. Such interest is motivated not only by poten- 
tial technical applications, but also on purely scientific 
grounds. Research in this field may lead to a deeper 
understanding of the fundamental working principles be- 
hind Nature's behavior, and often original new physical 
theories emerge. One of them, recently proposed in Ref. 
[3, predicts the existence of a class of systems whose crit- 
ical behavior lies outside the scope of the Landau theory 
of phase transitions^. Critical points in these systems are 
characterized by the deconfinement of fractionalized ex- 
citations, parameterizing the original degrees of freedom, 
which occurs right at the transition. It was observed 
that this scenario can, in principle, be realized in spin 
systems, which exhibit a second-order phase transition 
point characterized by the simultaneous breakdown of 
a continuous (e.g. spin SU{2)) and a discrete (e.g. lat- 
tice) symmetries, in such a way that symmetry groups on 
opposite sides of the transition are not group-subgroup 
related. Such critical points cannot be described in the 
framework of Landau's theory. 

According to Ref. [3, there should exist a substan- 
tial number of spin systems, which exhibit deconfined 
critical points. For instance, frustrated two-dimensional 
(2D) antiferromagnets (AFs), like the J1-J2 model, are 
believed to fall into this category. However, there seems 
to be no experimental or theoretical proof of this claim. 
Another class of models believed to display such a behav- 
ior includes non-frustrated AFs with multi-spin exchange 
interactions. One such model was studied by Sandviki 
and other authors^i^ and, although seemingly artificial, 
it provides a playground for testing new theories. Their 
Quantum Monte Carlo (QMC) simulations claimed nu- 
merical evidence for the deconfined quantum criticality 
scenario. 

The model studied in Ref. is a Heisenberg AF with 
a subset of four-spin ring exchange interactions, defined 



on a square lattice (named the J-Q model) 



(ijki) 

where Q > 0, denote sites in a 2D square lattice 

and Si are spin-1/2 operators. The first summation ex- 
tends over bonds (nearest neighbor sites). The second 
term contains two sums over plaquettes (sites of the dual 
lattice): first, (ij) and {kl) denote parallel horizontal 
links of the plaquette, and then {ik) and (jl) correspond 
to parallel vertical bonds. It was concluded^ that there 
exists a critical point at Qc/ J ~ 25 separating the an- 
tiferromagnetic phase from a valence-bond solid (VBS) 
state, whose nature is, strictly speaking, unclear— but 
the calculations suggested a columnar (dimer) order in 
this paramagnetic region. 

In the present paper we study the phase diagram of 
the J-Q model, using a recently proposed hierarchical 
mean-field (HMF) techniqueii^. The main idea of the 
method revolves around the concept of a relevant degree 
of freedom (a "quark") - spin cluster in this particular 
case - which can be used to build up the system. The 
initial Hamiltonian is then rewritten in terms of these 
coarse-grained variables and a mean-field approximation 
is applied to determine properties of the system. Thus, 
the (generally) exponentially hard problem of determin- 
ing the ground state of the system is reduced to a polyno- 
mially complex one. At the same time, essential quantum 
correlations, which drive the physics of the problem, are 
captured by the local representation. In other words, pro- 
vided the quark is chosen properly, even a simple single 
mean-field approximation, performed on these degrees of 
freedom, will yield the correct phase diagram. Moreover, 
our HMF ansatz provides an educated nodal surface that 
can, in principle, be used in conjunction with fixed-node 
(or constrained patb^) QMC approaches to further im- 
prove correlations, and thus energy estimates, in those 
cases when there is a sign (phase) problem. 

It is important to emphasize the simplicity of our 
method. In this work we concentrate on symmetries of 
the various phases, exhibited by the J-Q model. By using 
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Figure 1; Schematic phase diagram of the J-Q model, ob- 
tained in the present paper. The translationally invariant AF 
phase with broken SU{2) symmetry and the singlet param- 
agnetic state are separated by a quantum phase transition at 
Q = Qc- Results of our calculations indicate that this tran- 
sition is continuous, but we can not rigorously discard the 
possibility that it becomes weakly first order in the thermo- 
dynamic limit. The non-magnetic phase breaks the lattice 
translational symmetry and is a plaquette paramagnet. 



a more sophisticated variational ansatz (e.g. a Jastrow- 
type correlated wavefunction) , one can also improve nu- 
merical values of the observable quantities and phase 
transition points, but the physical picture will remain 
intact. Nevertheless, the HMF method was quite accu- 
rate to yield the quantitatively correct phase diagram^ of 
the J1-J2 model, whose behavior is driven by the inter- 
play of two gapless phases: the Nccl and columnar AF 
states. In the J-Q model the large-Q phase is gapped. 
Due to its real-space nature, the HMF method should be 
appropriate for this model. Indeed, our recent studiea^^ 
of another gapped system - the Heisenberg model on the 
Shastry-Sutherland lattice - support this assumption. 

Our findings are summarized in Fig. [TJ There indeed 
exists a non-magnetic phase, which we found to be of a 
correlated plaquette type, and not of a dimer character, 
separated from the Neel state by a second-order phase 
transition. However, the numerical value of the critical 
point, Qc/J ~ 2, which we obtained, is quite different 
from that of Ref. [J. Our results are consistent with 
data obtained from exact diagonalization of finite spin 
clusters which, given the fact that the system is gapped 
in the paramagnetic phase for Q ^ Qc, should be reliable 
in this region. Although we found the phase transition to 
be of the Landau second-order type, due to the real-space 
nature of the method (which explicitly breaks the lattice 
translational invariance), we cannot rigorously rule out 
the possibility that this phase transition becomes weakly 
first order as one scales the degree of freedom towards 
the thermodynamic limit. 

In the next section we set up the formalism. The re- 



sults are presented in Sec. 
in Sec. [iVl 



nil followed by a discussion 



II. COARSE GRAINING AND HMF 
APPROXIMATION 

For our purposes it is convenient to separate the two 
and four spin terms in the J-Q Hamiltonian: 



H = - 



2NQ 
16 



J+^)J2s,Sj 



(ij) 

Q {SiSj){SkSi). 
(ijki) 



(1) 



A satisfactory coarse graining procedure should parti- 
tion the lattice into spin clusters (quarks), containing Nq 
sites, that explicitly preserve symmetries of the Hamil- 
tonian. In particular, the J-Q Hamiltonian is explic- 
itly spin-S'f7(2) invariant. Moreover, it is invariant un- 
der transformations from the lattice rotational group C4. 
Therefore, we will consider only symmetry preserving de- 
grees of freedom: (i) plaquettes (2x2 spin clusters) and 
(ii) 4x4 spin clusters. Each cluster state will be asso- 
ciated with a hard-core bosonic operator 7. These oper- 
ators are Schwinger bosons of iS'C/(2^9) and must obey 

the local constraint: J^a'^ia'^ia = 1- They define the 
hierarchical language^ for our problem. 

From the form of Eq. ([T]) it is clear— that the (equiv- 
alent and exact) bosonic Hamiltonian will contain not 
only two-body scattering processes, but also four-boson 
interactions. Therefore, we can write down symbolically: 



i 

E( tt'2 \^'\^2 t t 

4 
4 

+ V ('ii-.4«r'i"3;4< TT t 

{iii2i3i4)v M=l 



(2) 



where a, . . . label states in the Hilbert space of a quark, 
i, j, . . . denote sites in the coarse grained lattice, and sum- 
mations are assumed over all repeated indices. The term 
with H?^^ encodes two-body interactions, while the last 
two lines describe the correlated four-boson scattering. 
The superscript h indicates that 1112 and 13*4 are hori- 
zontal links of a plaquette, and similarly v denotes the 
case when iii^ and ^2*4 are vertical links of the same 
plaquette. 

We will investigate the phase diagram of the J-Q model 
using the HMF approximation whose variational state as- 
sumes that the hard-core bosons form an insulating state. 
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Further, we introduce a new set of bosonic operators, re- 
lated to the old ones by a real site-independent canonical 
transformation: 



and write the variational ground state in the form: 

l'/'o) = nrlolo), 



(3) 



where n = denotes the lowest energy single-particle 
mode (we shall also denote: i?° = Ra), and |0) represents 
the vacuum. It is important to emphasize that although 
the coarse graining procedure preserves the symmetries of 
the Hamiltonian, some of them can be spontaneously bro- 
ken at the mean-field level as a result of self-consistency. 
In particular, the columnar dimer state is contained in 
the wavefunction ^ although, as we will see below, it 
never appears as a stable solution. 

We have explicitly separated the four-boson interac- 
tion in the Hamiltonian Q into horizontal and vertical 
link contributions. This distinction is important because 
these two terms must be properly symmetrized to fulfill 
bosonic statistics. In particular, the term H^^^ has to 
be symmetrized only with respect to indices in the same 
group, and groups as a whole (groups are separated by 
semicolons), i.e. one needs to take into account only the 
following permutations: (1 <-> 2), (3 <-» 4) and simultane- 
ously (1 ^ 3, 2 ^ 4). Analogously, in the term iJ^" only 
the permutations (1^3), (2 <-> 4) and (1 ^ 2, 3 ^ 4) 
should be accounted for. 

The problem then reduces to minimization of the en- 
ergy functional: 



N 



-{Hci)a'aRa'Ra + 

+(^i^t):::fii?<i?a.+ 

1^=1 

, , , , i 



(4) 



4h 
int 



under the constraint RaRa = 1, which leads to the self- 
consistent eigenvalue equation: 



HHF)^^Rb — ^J'Ra 



(5a) 



with the chemical potential /.i being the lowest eigenvalue 
of the Hartrce-Fock Hamiltonian: 

Once the amplitude Ra is determined, the ground state 
energy (GSE) can be computed using Eq. 
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Figure 2: The plaquette lattice. Thick lines denote interac- 
tions J -\-Q/2. The circles indicate four-spin terms of strength 
Q in the Hamiltonian Small-sized numbers label spins 

within a plaquette, while the larger ones label plaquettes. 



Although we have formulated the HMF method for 
spin-1/2 systems, it can be straightforwardly extended 
to higher spins as well. 

Besides the GSE we will also be interested in com- 
puting the staggered magnetization M^, and the two- 
component VBS "order parameter"—: 



X 



which allows us to characterize lattice point group sym- 
metries of a state. 

In the rest of this section we will sketch the HMF cal- 
culation of E'o; and ^ for the case of plaquettes, and 
only present final expressions for the 4x4 clusters. The 
interested reader is referred to Ref. 0, where the tech- 
nique is analyzed and developed in greater detail. For 
simplicity we shall put J = 1. 



(i) The plaquette degree of freedom 

We start by considering the simplest way to cover the 
lattice - with plaquettes, as shown in Fig. [51 At the same 
time we introduce notations and concepts, which will be 
used in the following subsection. 
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Figure 3: Connectivity of the lattice formed by 4 x 4 spin clus- 
ters. Small circles indicate spins. Bold numbers label the 4x4 
clusters, while thin ones denote plaquettes. Other notations, 
are the same as in Fig. [2] 



The Hamiltonian for an isolated plaqucttc has the 
form: 

Ha = (^1 + 1^ (5l + S4) (52 + ^3) - 

~ Q[{SiS2){S3S4) + {SiSs){S2Si)]. (6) 

The interaction of this plaquette with the rest of the sys- 
tem can be conveniently partitioned according to ^ as: 

rr tt2 , Tjih , Tjiv 

-Hint — ^int + ^int + ^int: 

where appropriately symmetrized individual terms are 
given by: 

HL = [(^11 + Su){S22 + ^23) + 

+ (5l2 + Sl3)(521+524)]- (7a) 

[(>S^ll5'l2)(S'23S'24) + (S^llS'23)(S'l2S'24) + 

+ (S'i2S'l4)('S'21<S'23) + («S'l2'S'2l)(S'l4S'23) + 
+ ('S'l35'l4)(S'2l-S'22) + (S'i3S'2l)(S'i4S'22) + 
+ {SllS22){S 13'S'24) + (S'll'S'l3)(S'22S'24)] ; 

^^hit ~ [(5^145^23 + 'S'24S'i3)(5'32S'4i + 531542) + 

o 

-|-(534543 + 544533)(5i252i + 5ii522)] ; (7b) 



H^nt — [('S'14532 + 512^34) (523 541 + 52l543) + 

-|-(524542 + 522544)(5i353i + 5ii533)] . (7c) 

It is convenient to work in the basis which diagonal- 
izes the Q- independent part of iJp, Eq. ([5]). Such is, 
for instance, the basis of eigcnstates of the total angular 
momentum of the plaquette: 

\a) = \hl2LM), (8) 
= 5i + 54; ^2 = 52 + 53; L ~ li + l2- 

The matrix elements, which appear in Eq. ^ 
{Hl^';^ Ma[a',\Hl,\a,a2); 

(^int )aia.;a3a4 = <"l°2;a3a4l^^int \aia2; a3a4) 

can now be computed using the angular momentum ad- 
dition theorems. 

The staggered magnetization (along the 2;-axis) within 
a plaquette is given by 

= + 5| - S^, - SD^^RaRb, (9) 

while the function 5* can be written in the plaquette rep- 
resentation as: 

Re* =^ ^ [5,1 5,2 + 5,35,4] - 

i 

~ [Si2Si^x^i + Si4Si^x,3] ; (10) 

i 

Im* =^ E + '^'2S'.4] - 

i 

— ^ [5,;35i+y,i + 5i45,+y^2] ■ 

i 

In these equations the indices i and x denote sites and 
basis vectors of the plaquette lattice. 

(a) 4x4 spin clusters 

The coarse grained lattice obtained by choosing the 
4x4 cluster as a degree of freedom is shown in Fig. [S] 
Each spin operator carries three indices: label of a clus- 
ter, label of a plaquette within this cluster, and the posi- 
tion within this plaquette. Writing down the cluster self- 
energy and the inter-cluster interactions is a straightfor- 
ward, but tedious task, which can be accomplished along 
the lines presented in the previous subsection. Therefore, 
here we give only final expressions for AI^ and 'f. 

The staggered magnetization of a cluster is given by 
an equation analogous to ([5]): 

4 

M, = ^(5fi + - - S:r,)^,^RA>RA, (11) 

i=l 
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Figure 4: Ground state energy as a function of Q/J for pla- 
quettes [case (i)] and 4x4 spin clusters [case (ii)]. 



where the summation extends over plaquettes within a 
cluster, and we used capital indices A to label states of 
the cluster. The function 'I' can be written as 

R.e^' = — '^^{SiiiSii2 + Sii3Sii4 + Si3iSi32 + 

i 

+ SiS^SiSi + Si2lSi22 + Si23Si24 + 
+ Si4iSu2 + Si43Si44 — Sii2Si21^ 
— Sii4Si23 — Si32Si4i — Si34Si43) — 

— ^ [Si22Si+x,ll + Si24Si+x^i3 + 

i 

+ Si42Si+x^3i + Si44Si+x^33\ ; (12) 

Im* = ^ ^(S'jiiS'jia + Sa2S^i4 + SaiSi- 



i23" 



+ Si22Si24 + Si3iSi33 + Si32Si34 + 
+ Si4iSi43 + Si42Si44 — SnaSisi — 
— Sil4Si32 ~ S 1238141 — 81248142) — 

^ [Sa^Si + y^ll + 8iZ4Si + yS2 + 



+ 8i43Si+y^21 + Si448iJ^y^22\ ■ 



III. RESULTS 

We can now proceed with solution of the mean-field 
equation (HJ), supplemented by Eqs. Q for the case 
[i) and analogous expressions in the case (ii). 

The physical quantities that we want to compute in the 
first place are the GSE and the staggered magnetization. 
These are given by Eq. ([4]) and Eqs. ([9]) for the case (i), 
and (fTT|) for the case (ii). In Fig. [3]we present GSEs for 
both degrees of freedom. All energies monotonically de- 
crease with increasing Q/ J as a consequence of the nega- 



Figure 5: Second-order derivative (PEo/NdQ^ (main panel) 
as a function of Q/J for cases (i) and (ii). The discontinuity 
at Q/J ~ 1.61 (2 X 2) and Q/J - 2.0 (4 x 4) indicates a 
second-order phase transition point. The inset shows the ex- 
trapolation of the jump g — J / N dQ'^\'Q''^ to Nq 00. 



tive sign in front of the last term in Eq. ([T]). At some crit- 
ical value of Q = Qc the system undergoes a phase tran- 
sition from the Neel state at small Q to a spin-disordered 
state at Q > Qc- This transition can be seen either from 
the second derivative of the GSE, (PE^/NdQ^ , shown in 
Fig. [5l or from the staggered magnetization as a func- 
tion of Q/J, presented in Fig. [6l Using these plots one 
obtains the numerical values Qd J = 1-61 for plaquettes 
and Qd J = 2.00 for 4x4 clusters. Although the jump 
g = JcPEt:,/NdQ'^\'Q'd^Q is numerically small, it remains 
finite: g — > 0.016, if extrapolated to the thermodynamic 
limit, based on these two points (see the inset to Fig. [5]). 




Figure 6: Staggered magnetization (main panel), Eqs. (|9]) 
and pif) . as a function of Q/J for cases (i) and (ii). The 
values of Qc are: Qd J = 1.61 for case (i) and Qd J = 2.00 
for case (ii). The inset shows the scaling of Qc- 
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Figure 7: The "order parameter" for the two cases, studied 
in this paper. Notice the coincidence of curves for Re*]/ and 
Im*!'. For Q ^ Qc tliis implies the plaquette nature of the 
quantum paramagnetic state. 

The finite-size scaling of the critical point itself, presented 
in the inset to Fig. [SI shows that hniAr_j_oo Qc/J = 2.13. 
In order to demonstrate that our results are reliable, 
wc compute limiting values of the GSE and the mag- 
netization at Q = 0: limAr^_too -E-o/^^^ = —0.64 and 
limAr^_,tx) Mz = 0.39. These numbers should be com- 
pared to the accepted QMC resultsi^: Eq/NJ = -0.67 
and Mz = 0.31. We note, finally, that due to few data 
points, the finite-size scalings presented here are quali- 
tative, and are intended to provide only an estimate for 
the extrapolated quantities in the thermodynamic limit. 

Let us now discuss the symmetries of the various 
phases. The antiferromagnetic state, which occurs for 
Q < Qc, is known to preserve the lattice rotational sym- 
metry C4, and spontaneously breaks the spin SU (2) sym- 
metry. The nature of the paramagnetic phase, stabilized 
for Q ^ Qc, can be unveiled by computing expectation 
values of the function 4" given by Eqs. PT7)1 and p2p for 
cases (?) and (ii), respectively. Although 5* is an integral 
quantity, it is sufficient for the purpose of discriminat- 
ing between plaquettized and dimerized ground states. 
Namely, a plaquette phase preserves the four-fold lattice 
rotational symmetry, implying 

Rc* = Im*, (13) 

while in a dimerized state this equality does not hold. 
In Fig. [3 wc show Re* and Im*. The equality (fH)) 
is satisfied throughout the phase diagram. This fact is 
not surprising in the antiferromagnetic phase, but in the 
paramagnetic region it presents a strong evidence against 
any type of dimerized ground states. Although such 
states were allowed in the process of minimization, the 
C4-symmetric states always had lower energy. In fact, 
the ground state in the non-magnetic region is a plaque- 
tte paramagnet, with each plaquette being in its singlet 



Figure 8: Contributions to the GSE from the J and Q terms in 
([T]) for a 4 X 4 spin cluster with periodic boundary conditions. 
The unimportant term —NQ/8 is omitted. 



ground state. However, due to the tensor nature of in- 
teractions in ^ these plaquettes are interacting. 

As already mentioned in the Introduction, the coarse 
graining procedure explicitly breaks the lattice transla- 
tion invariancc, which should be restored in the thermo- 
dynamic limit. Extrapolation to A^^ — > 00 shows that 
Re^fjImVlf — > —0.04, suggesting that the translation in- 
variance is indeed being recovered. 



IV. DISCUSSION 

Our calculations, presented in the previous section, 
demonstrate that the Hamiltonian ([T]) exhibits a phase 
transition point separating the Neel ordered state from a 
paramagnetic phase with broken translational invariance, 
in agreement with conclusions of previous works^i^. 
Most importantly, besides establishing the existence of 
a phase transition, we were also able to unveil the nature 
of the paramagnetic phase and show that a correlated 
plaquette state is favored over a columnar dimer state 
which, although not conclusively, seems to be preferred 
according to previous calculations^. 

However, despite qualitative agreement, there is a 
quantitative discrepancy in the numerical value of Qc- 
Namely, the value obtained in the present paper is much 
smaller than the one presented in Ref. 0. Although we 
cannot provide a rigorous explanation for this discrep- 
ancy, we would like to make some qualitative remarks in 
the following: 

First of all, it is clear that the variational wavefunc- 
tion ([3]), being a low-density ansatz, generally leads to 
an under-estimation of the four-boson scattering terms 
in Eq. In order to understand how significant this 

error is and check that the results, presented in Figs. [3HS1 
are reasonable, wc used data from the exact diagonaliza- 
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tion of 4 X 4 spin clusters to compare magnitudes of the 
two terms: the ones, proportional to J and Q, in Eq. ([1]). 
On physical grounds one would expect a phase transition 
to occur when these terms become comparable. Figure [5] 
presents the two contributions and their dependence on 
Q/J. Of course, the crossing point at Q/J ^ 1 does not 
determine the critical value Qc, but it provides a clue on 
where the phase transition may occur. Since the system 
is gapped in the paramagnetic phase, one can argue that 
the size 4 x 4 is large enough to describe the thermody- 
namic limit. Indeed, QMC data for Q/J = 10 indicates 
that the GSE converges very rapidly with increasing sys- 
tem sizai^. Also, calculations analogous to that shown 
in Fig. m performedi^ for systems up to 16 x 16 sites, 
indicate that the magnitude of the crossing point stays 
of order unity. 

Second, we would like to emphasize that, although 
there is no question about the correctness of the QMC 
studies of Refs. [HsHolfl^. the procedure used to extract 
physical quantities, like Qc, from the raw statistical data, 
is not straighforward and requires certain assumptions^. 
Therefore, it is desirable to have another independent 
determination of the phase transition point, for example, 
from the data on staggered magnetization, computed in 
the entire Q range. While such calculations would def- 
initely help to resolve this issue, surprisingly, they have 
never been performed. QMC computations of finite lat- 
tices does not suffer from the infamous sign problem in 
this case, thus it yield better energy and magnetization 
values than the ones obtained here. Our approach, on the 
other hand, focuses on establishing symmetry properties 
of different phases, rather than improving numerical val- 
ues for observable quantities. It is this fact, which en- 
ables us to detect phase transition points within a simple 
framework. 



Our conclusions raise another important question re- 
garding the nature of the phase transition. We find it to 
be of the Landau type. Although the finite-size scaling 
of the second-order derivative of the GSE, presented in 
the previous section, displays a finite jump as Nq — > oo, 
there is no way to rigorously prove it. Thus, the pos- 
sibility of a weakly first order transition at Qc cannot 
be completely excluded. Indeed, in Ref. [ij it was ar- 
gued that this phase transition, which was claimed to 
occur at the same point as in Ref. 0, is of the first or- 
der. As any real-space method our approach explicitly 
breaks translational invariance, and although the finite- 
size scaling for 4* implies that this property is restored 
with increasing cluster size, we cannot provide a rigorous 
symmetry-based analysis. 

In summary, we determined the phase diagram of the 
J-Q model ([1]), by using the recently proposed hierarchi- 
cal mean-field approach^. It was shown that there exists 
a single (i.e. universal) mean-field framework (variational 
ansatz for the ground state), which gives the complete 
phase diagram of the model. In particular, we found that 
there exists a critical point at Qc ~ 2. 13 J, which sepa- 
rates the antiferromagnetic phase from the non-magnetic 
state. The latter breaks lattice translational invariance 
and was shown to represent a correlated plaquette para- 
magnetic phase. Our results suggest that the phase tran- 
sition at Qc is of a Landau second order type, even in 
the thermodynamic limit Nq —^ oo, although we cannot 
rigorously exclude the possibility for it to become weakly 
first order. 

We are indebted to A. W. Sandvik for bringing this 
problem to our attention and numerous discussions and 
private communications in the course of completion of 
this work. JD acknowledges support from the Spanish 
DGI under the grant FIS2006-12783-C03-01. 
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